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1— ( ' Abstract. We review recent developments in the theory of interacting quan- 

^ , tum many-particle systems that are not in equilibrium. We focus mainly on the 

\L{ • nonequilibrium generalizations of the flow equation approach and of dynamical 

Ti ' mean- field theory (DMFT). In the nonequilibrium flow equation approach one 

. ; , first diagonalizes the Hamiltonian iteratively, performs the time evolution in this 

C^ ■ diagonal basis, and then transforms back to the original basis, thereby avoiding a 

C ' direct perturbation expansion with errors that grow linearly in time. In nonequi- 

I , librium DMFT, on the other hand, the Hubbard model can be mapped onto a 

'■^ ■ time-dependent self-consistent single-site problem. We discuss results from the 

^ ' fiow equation approach for nonlinear transport in the Kondo model, and further 

O . applications of this method to the relaxation behavior in the ferromagnetic Kondo 

,^ , • model and the Hubbard model after an interaction quench. For the interaction 

quench in the Hubbard model, we have also obtained numerical DMFT results 

I , using quantum Monte Carlo simulations. In agreement with the flow equation ap- 

^ ■ proach they show that for weak coupling the system relaxes to a "prethermalized" 

[*--. ' intermediate state instead of rapid thermalization. We discuss the description of 

^\ , nonthermal steady states with generalized Gibbs ensembles. 

O' 
vn ; 

jy-^ ! 1 Introduction 

o; 

^^ . Strongly correlated electron systems and their rich phase diagrams continue to play a central 

^^ ' role in modern condensed matter physics. Key theoretical developments were the solution of 

^ , the Kondo model as the paradigm for correlated quantum impurity models using the numer- 

V"! ' ical renormalization group J^, and the solution of the Hubbard model as the paradigm for 

rS I translation-invariant correlated electron systems within dynamical mean-field theory (DMFT) 

3 . [2l3j . Until about ten years ago these investigations nearly exclusively focussed on equilibrium 

■ ■ ■ ' or near-to-equilibrium (linear response) situations. This was mainly due to the fact that ex- 

periments probed either equilibrium or linear response properties. For example electrical fields 
applied to bulk materials are typically too weak to drive a system beyond the linear response 
regime. 

The experimental situation changed completely in the past ten years due to three key 
developments. First of all, the realization of various model Hamiltonians using cold atomic 
gases paved the way to study the real-time evolution of quantum many-body systems and 
quantum quenches away from equilibrium. The seminal experiment in this context was the 



demonstration of collapse and revival oscillations in an ultracold gas of rubidium atoms that 
are described by a Bose-Hubbard model quenched from the superfluid to the Mott phase by 
Greiner et al. [3]. The second modern experimental development is femtosecond spectroscopy 
[5I6I7I8) , which permits to study the electron relaxation dynamics in pump-probe experiments. 
Third is the realization of correlated quantum impurities in Coulomb blockade quantum dots 
[9I10I11J . which can easily be driven beyond the linear response regime with moderate applied 
voltage bias [12 . 

These experimental developments have led to numerous theoretical investigations of 
nonequilibrium quantum many-body systems in the past decade. Still it is fair to say that the 
nonequilibrium properties of correlated systems are much less understood than their equilibrium 
counterparts. The main reason for this is the lack of reliable theoretical methods that can cope 
with the double challenge of strong correlations and nonequilibrium situations. In this paper 
we therefore highlight two such theoretical approaches that have been applied successfully to 
nonequilibrium quantum many-body systems, namely the analytical flow equation method [13) . 
generalized to nonequilibrium IT , and nonequilibrium dynamical mean- field theory [151161 . We 
discuss how these methods contributed to a first understanding of some paradigms of correlated 
electron physics in nonequilibrium. 

The paper is organized as follows. In Sec. [2] we introduce the fiow equation approach for an- 
alytic diagonalization of quantum many-body systems, especially its generalization to nonequi- 
librium real-time evolution problems. Sec. [3] deals with the application of this approach to the 
ferromagnetic Kondo model, which can thereby be solved in a controlled way even for asymp- 
totically large times. The results are in very good agreement with numerical methods [T7] and 
also establish a key mechanism for thermalization after an interaction quench in the Hubbard 
model. The time evolution of the Hubbard model is then analyzed in Sec.|4]using the fiow equa- 
tion method. Sec. [5] contains an application of the flow equation method to a second important 
class of nonequilibrium problems, namely to transport beyond the linear response regime. In 
Sec.|6]wc briefly review nonequilibrium DMFT, which maps the lattice system onto an effective 
single-site problem with a self-consistency condition. We show how this self-consistency condi- 
tion can be reduced to a single equation in some cases, which can reduce the numerical effort 
for the solution of DMFT, in particular in the absence of translational invariance in time. In 
Sec. [7] we discuss DMFT results for interaction quenches in the Hubbard model, which agree 
very well with the results from the flow equation method for small values of the Hubbard inter- 
action (Sec. S]). For quenches to intermediate Hubbard interaction the system thermalizes on 
short time scales, which shows that correlated systems in isolation can reach a new equilibrium 
state, not because of a coupling to external baths but due to the interactions between particles. 
In Sec. [5] we compare this behavior to that in a special solvable Hubbard model, for which 
the system tends to a nonthermal state due to its integrability, i.e., the conservation of many 
constants of motion. Finally in Sec. |9] we discuss the concept of generalized Gibbs ensembles, 
which can make statistical predictions both for integrable and nearly integrable Hamiltonians. 

2 Flow equation approach to quantum time evolution 

Stationary eigenstates of a many-particle Hamiltonian are fundamental for the discussion of 
quantum many-particle systems in equilibrium. A typical class of nonequilibrium situations are 
systems that are prepared in a quantum state \^i) which is not an eigenstate of the Hamiltonian 
H that drives its time evolution. In this case, observables that do not commute with H will 
generally become time-dependent. The canonical way to evaluate the real-time evolution of such 
observables is the Keldysh technique. One of the notorious difficulties with this approach is that 
often the limit of weak interactions does not commute with the limit of long times. A partial 
sum over certain diagrams is usually not sufficient to guarantee a controlled approximation in 
the limit of long times. 

Another approach for calculating time-dependent observables is Heisenberg's equation of 
motion for an observable A, 
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Fig. 1. The forward-backward transformation scheme induces a non-perturbative solution to the 
Heisenberg equation of motion for an operator A. U denotes the fuU unitary transformation that 
relates the B = to the B = oo basis 1201. 



Again, a direct perturbative expansion of this equation of motion, e.g., in the interaction 
strength, is usually not controlled in the limit of long times. However, one might think of 
changing the basis representation of Eq. ([T]) . Ideally, this basis representation would transform 
H into a non-interacting form. This idea is precisely at the heart of the flow equation approach, 
invented by Wegner in 1994 [T^, and independently by Glazek and Wilson [18119) . 

The main concept is to implement this transformation as a sequence of infinitesimal unitary 
transformations that allows to separate energy scales during the process of transformation. 
This is achieved by introducing an appropriate generator r](B) that parameterizes a family of 
unitarily equivalent Hamiltonians. By solving the differential equation 






(2) 



with the initial condition H{B = 0) = H , a unitary equivalent family of Hamiltonians is con- 
structed. Wegner's ingenious choice of a canonical generator is to construct it as the commutator 
of the diagonal (i?o) with the non-diagonal part (ifint) of the Hamiltonian, 
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This ensures that H{B) becomes more and more energy-diagonal with an effective band width 
ylfeq during the flow. One can easily verify the identification ylfoq ex B~^^'^. Generically one 
reaches a non-interacting diagonal form in the limit B ^ oo [13] . 

Intuitively, the basis defined by constructing the non-interacting normal form of the Hamil- 
tonian should be much more suitable to solve Heisenberg's equation of motion. This has been 
used in numerous examples for evaluating equilibrium correlation functions on all energy scales 
using flow equations |14j. At the same time it also allows the calculation of nonequilibrium 
real-time evolution problems. Indeed, the analogy of constructing normal forms of interacting 
Hamilton functions in order to integrate Hamilton's equation of motion is a very successful con- 
cept in classical mechanics, known as "canonical perturbation theory" [21,. Based on the flow 
equation approach, this idea can be directly applied to quantum many-body systems [20I22J . 
The general setup is described by the diagram in Fig. [I] where I'Fi) is some initial non-thermal 
state whose time evolution one is interested in. Here and in the following, we describe operators 
and coupling constants in the B = oo basis by a tilde. In order to study the real-time evolution 
of a given observable A that one is interested in, the observable is transformed into the diagonal 
basis by solving the differential equation 
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with the initial condition 0{B = 0) = A. The key observation is that one can now solve the real- 
time evolution with respect to the energy-diagonal H exactly, thereby avoiding any errors that 
grow proportional to time (i.e. secular terms): this yields A{t). Now since the initial quantum 



state is given in the B — basis, one undoes the basis change by integrating Eq. ^ from B — oo 
to B — (backward transformation) with the initial condition 0{B — oo) = A{t). One therefore 
effectively generates a new non-perturbative scheme for solving the Heisenberg equation of 
motion for an operator, A{t) = e*^*A(0)e~*^*, in exact analogy to canonical perturbation 
theory. In a first successful application, this approach has been applied to dissipative quantum 
systems [22|20J . In this paper, we will discuss recent work on the ferromagnetic Kondo model 
(Sec. [3]), the fermionic Hubbard model (Sec. H]) and on steady state transport through Kondo 
dots (SecEl). 

Usually, the implementation of the flow equation approach relies on approximations that 
allow to truncate the hierarchy of flowing interaction terms generated during the transformation 
of an interacting many-particle Hamiltonian and all observables. Systems that are accessible by 
perturbative RG are ideally suited for that purpose, since they allow for a controlled expansion 
around a weak-coupling fixed point. 



3 Ferromagnetic Kondo model 

A well-known example for a weak-coupling many-particle system can be realized in the ferro- 
magnetic regime of the Kondo model 

^ = E '"^L^ka + E 4fc^^4'fe + E Jt'kis+s^,^ + s-s+j , (5) 

ka k,k' k,k' 

where S = S^ztiS^ and likewise for the conduction electron spin densities. In the following, we 
consider constant ferromagnetic exchange couplings with \J'^\ < J". The canonical generator 
is immediately obtained from ri{B) = [Ho,Hint{B)], where the flowing interaction Hint{B) 

is parametrized by the flowing couplings ^^^/(-B) and J^k'^^)- ^^ ^^^ Fermi surface (with 
Jj^ '"^ = J^'ll and the density of states p), it can be shown that these couplings reproduce the 
usual poor man's scaling equations [23] 
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For ferromagnetic couplings, the non-interacting fixed point of Eq. (j5]) is stable and allows 
for a controlled perturbative expansion of all flow equations used to transform H and S^. 
This example allows to analytically implement the approach outlined in Sec. [2] to calculate the 
time-dependent impurity magnetization {S^{t)) with respect to the initial state 

1^,) = I t) ® |FS) . (7) 

The state \'Fi) has the physical interpretation of an initially decoupled system of a polarized 
impurity spin | f) and a non- interacting Fermi sea in equilibrium, |FS). The remaining technical 
steps consist of transforming the operator S^ according to Fig. [U in order to evaluate the 
observable {S%t)). 

The impurity magnetization is obtained as follows. It is straightforward to work out the 
ansatz for the flowing spin operator as 

S'-{B)^h{B)S^+Y,Jk'k{B):{Sxsk'ky: • (8) 

kk' 

In the limit B — )■ oo, the coupling constant h{B) approaches the value h{B -^ oo) = 1 + pJ/2 + 
0{J^) if j" = J"*". The value h{B — > oo)/2 is equal to the equilibrium impurity magnetization 
in presence of an infinitesimal Zeeman term 0+5^ and has been calculated to the same accuracy 



by Abrikosov 24 by a summation over parquet diagrams. According to our strategy for the 
solution of Heisenberg's equation of motion, the nonequihbrium magnetization follows from an 
ansatz 

S'{B, t) =. h{B, t)S' + Y, jk'kiB, t):{Sx s^kY ■ ■ (9) 

fcfc' 
Obviously, the flow equations for this operator have the same form as the time-independent flow 
equations for the ansatz in Eq. ([5]). It is furthermore trivial to determine the initial condition 
as S^(B ^oo,t) = h{B -^ oo)S' + Y.kM' lk'k{B -^ oo)e'*(^'='-^'=) : {S+s-,^ + ^"4 J : . Up 
to neglected normal ordered terms of 0{J^), the operator S^{t) readily follows from integrating 
Eq. ^ with the initial condition S^{B — > oo,t). Using these steps, the formal result for the 
impurity magnetization reads 

iS'it)) = ^ + E % (e'*(^'=-^'='^ - l)nik')[l n{k)] . (10) 

By solving the flow equations for the couplings ^kk> for momenta close to the Fermi surface, it 
is therefore possible to directly obtain the long-time dynamics of the impurity magnetization. 
More details of the calculation can be found in Refs. |25I17| . For isotropic couplings and using 
a flat band with the dimensionless range [—1, 1], this yields the asymptotic behavior 

<^' W> ^ I (r^F^ + 1 + pJ + o{j')^ . (11) 

This asymptotic behavior is confirmed by numerical calculations depicted in Fig. [21 For 
anisotropic couplings, the calculation is completely analogous and yields the result 



with the dimensionless coupling j" = pVjII^ — J^2 g^jj,^ q, r^ pj\\ j-qj- |j-L^j||| > 2. These 
results allow for two interesting observations: (i) the time-dependent contribution to {S^{t)) 
is described by jl'(yl = 1/t), where J"(^) fiows according to the scaling equations, (ii) The 
asymptotic result {S^{t -^ oo)),j/; — 1/2 = 2{{S^)eq — 1/2) is obeyed. Therefore, a factor of two 
distinguishes the asymptotic results for the equilibrium magnetization vs. the nonequihbrium 
magnetization. 

In this section we discussed a rare example where the nonequihbrium dynamics of an inter- 
acting many-particle system can be described analytically. The versatility of the flow equation 
approach to describe both equilibrium and nonequihbrium quantities with common scaling 
transformations allows to directly relate nonequihbrium relaxation laws to equilibrium prop- 
erties obtained by a conventional renormalization group calculation. In our example, we could 
establish that nonequihbrium dynamics is set by the fiowing coupling constants, evaluated at 
the energy scale equal to the inverse time scale of the relaxation process. Furthermore, we 
could show that an equilibrium and a corresponding nonequihbrium observable in a steady 
state differ by a factor of two. This observation appears in a much wider class of problems, 
e.g. in a class of discrete models, where this result can be stated in form of a theorem [26]. 
The same factor two also appears for the real-time evolution in the weak coupling phase of the 
quantum Sine-Gordon model |27j . where it connects the equilibrium with the nonequihbrium 
mode occupation numbers. It will be one of the central features for the real time evolution of 
the Hubbard model in Sec. [H since it defines the prethermalization and ensuing thermalization 
behavior. 



4 Small interaction quenches in the Hubbard model 

The forward-backward transformation scheme as it is depicted in Fig. [1] became a fruitful 
tool also for examining the nonequihbrium properties of closed translation-invariant interact- 
ing quantum many-body systems. Since the trace of the squared density operator p does not 
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Fig. 2. Time-dependent magnetization in the isotropic ferromagnetic Kondo model. D denotes the 
flow equation data points, which agree very well with time-dependent numerical renormalization group 
calculations depicted by x [T7]. The deviations for larger values of the coupling constant J can be 
understood from the perturbative nature of the flow equation calculation [23 • Using our analytical 
result we fitted the flow equation data against {S^{t)) = (1 -f aJ + [ln(i) — l/(cJ)]~^)/2 using a, c as 
fit parameters (lines). 



change under unitary time evolution, a closed quantum system which has been initialized in a 
pure state (defined by Tr[p^] = 1) will never relax to a thermal state (defined by Tr[p^] < 1). 
Nonetheless, the question of relaxation and thermalization can be addressed for expectation val- 
ues of observables. As the conventionally studied simple observables typically only probe one- or 
two-particle aspects of the many-particle quantum state, there is the generic expectation that 
the full unitary evolution cannot be resolved; then the apparent dynamics of these observables 
becomes arbitrary close to that of a thermal state. This is, however, not the case for integrable 
models where the dynamics is limited by additional conserved integrals of motion. In the fol- 
lowing we will show that the Hubbard model in more than one dimension, which is commonly 
considered nonintegrable, exhibits the phenomenon of prethermalization. It is characterized by 
the emergence of different relaxation times for various observables and has previously been 
observed numerically in cosmological models |28j . This implies that there exists a transient 
time regime during which the expectation values of some observables have already relaxed to 
their long-time values while others have not. The relaxation of the latter ones, however, is only 
delayed to a later time scale. 

We start with imposing nonequilibrium initial conditions by a quantum quench, i.e. a sudden 
change in a parameter of the Hamiltonian: We start in the ground state of a noninteracting Fermi 
gas and assume that the strength of the Hubbard interaction remains small such that we always 
remain within the Fermi liquid phase of the model. Therefore we first test the nonequilibrium 
properties of a Fermi liquid beyond Landau's theory. Quenches outside of the Fermi liquid phase 
are discussed using the DMFT approach in Sec. [T] Here we apply the forward-backward scheme 
which allows for an analytical ab initio real-time analysis of the subsequent relaxation dynamics 
of the kinetic energy and the momentum distribution |29l30j . It will turn out that the Fermi 
liquid picture provides a most suitable framework to understand the origin of prethermalization 
on the grounds of an analytical calculation. In the following paragraphs we will explain how 
prethermalization emerges from the interplay of the Pauli principle with quantum correlations 
induced by the two-particle interaction. 

The time evolution of a closed many-particle quantum system is fully described by its 
Hamiltonian. Let us assume that the latter contains different interactions between the various 
degrees of freedom on well-separated energy scales, for instance a two-particle interaction be- 



tween electrons, an additional electron-phonon coupling or one to an external light field. Then 
one trivially expects that the relaxation of a generic excited state passes through a sequence 
of transient time regimes: each of them corresponds to an energy scale of the Hamiltonian and 
exhibits its own characteristic signatures; only afterwards full thermalization of all degrees of 
freedom may be reached. 

Prethermalization of a Fermi liquid describes a similar sequence of time regimes which orig- 
inates, however, from a single two-particle interaction term in the Hamiltonian. Moreover, it is 
reminiscent of a common effective approach for weakly interacting systems which discusses the 
Hamiltonian time evolution in terms of two different relaxation mechanisms: dephasing of the 
initial state and two-particle scattering between unrenormalized momentum modes. This cor- 
responds to two distinct relaxation times: short time inelastic energy relaxation and long time 
elastic momentum relaxation. Prethermalization implies that these time scales separate and are 
observable in the relaxation behavior of different observables: 'bulk' quantities like the total ki- 
netic or potential energy of the system are 'momentum mode averaged' quantities which already 
relax due to dephasing. This can be seen as an implication of energy-time uncertainty which 
allows for rapid energy exchange at short times. However, in a translationally invariant system, 
(quasi-) momentum is conserved and the momentum distribution can only relax by momentum 
exchange due to scattering processes. Yet in a zero temperature Fermi liquid, scattering pro- 
cesses at the Fermi energy are suppressed by phase space restrictions which are imposed by the 
Pauli principle. Hence the momentum distribution represents the simplest example of a 'mo- 
mentum mode quantity' which relaxes on a later time scale. The separation of both time scales 
opens the transient regime of prethermalization during which a quasi-equilibrium description 
based on temperature cannot be defined for momentum mode quantities. 

We develop this scenario starting from a noninteracting Fermi gas at zero temperature, i.e. 
in the noninteracting ground state |^o)- Its momentum distribution exhibits a discontinuity 
of size Z° = 1 at the Fermi energy. Then the two-particle Hubbard interaction is switched on 
instantaneously, 

H = - Y. *^Aa^J- + UO{t){n,,^ - i)(n,; - i) , (13) 

where ci denotes a local fermionic destruction operator, rii = c\ci the number operator in real 
space (rifc in momentum space) and 0{t) the Heaviside step function. We assume that the 
nearest-neighbor hopping tij is larger than the local on-site two-particle interaction U < tij = 
1 to permit a weak-coupling description. The sudden change in the Hamiltonian initializes 
the system in a highly excited state and violates the fundamental prerequisite of Landau's 
equilibrium theory of a Fermi liquid, namely the adiabatic continuity of the noninteracting 
and the interacting Fermi system. Nonetheless we observe that important aspects of Landau's 
theory, in particular the quasiparticle picture of elementary excitations, are retained during the 
subsequent nonequilibrium dynamics. Applying the unitary perturbation scheme of Fig.[T]to the 
creation operator in momentum space implements its Heisenberg equation of motion. Thereby, 
it exhibits the decay of a physical fermion c^ into a superposition of various many-particle 
excitations which represent the 'dressing' of the noninteracting particle with particle-hole pairs 
due to interaction effects, 

Cka{t = 0, B) ™° Ckait, B) = hka{t, B) Cka + Mpqraaa{t, B) c],^cJg.Cr5- -|- . . . . (14) 

This ansatz is analog to the ansatz for the impurity magnetization ([5]) in the ferromagnetic 
Kondo model. It transfers the time evolution of operators (and/or their flow under infinites- 
imal unitary transformations) to a set of coupled differential flow equations for the time and 
i3-dependent prefactors hka{t, B) and Mpqrass{t,B). Higher order terms are truncated in a 
systematic way: only terms are kept which are relevant for a second order in U result of the 
momentum distribution function Nk{t) = [Dq] nk{t) |/2o)- At the Fermi surface, the coefficient 
hkpcrit) = \/ Z{t) relates to the time dependent value of the quasiparticle residue Z. The later 
mirrors the discontinuity of the momentum distribution function at the Fermi energy which 
is reduced under the time evolution. It approaches - in a formal long-time limit - a nonvan- 
ishing value ^^eq .^^Jj^jqJ^ rnismatches the corresponding equilibrium value Z^^^ by a factor 



/i = limt:^i/u{Z° - Z^^^{t))/{Z° - Z^^^) = 2. A comparison with exactly solvable models, 
e.g. the sudden squeezing of a harmonic oscillator |26', indicates that the mismatch is a nonper- 
turbative effect of nonequilibrium dynamics, while its numerical value /i = 2 is a perturbative 
result in the weak-coupling limit. Such a mismatch is the generic behavior of many systems and 
observables |26| . We have found it also comparing the nonequilibrium and equilibrium magne- 
tization in the ferromagnetic Kondo model (cf. Sec. [3]). For the Fermi liquid, the persistence of 
a nonvanishing quasiparticle residue up to late times has been confirmed recently [31]. Hence, 
up to a factor the momentum distribution still resembles that of a zero temperature Fermi 
liquid in equilibrium; therefore we conclude that a quasiparticle picture remains applicable. 
Note that, due to the mismatch, the corresponding quasiparticle momentum distribution is a 
nonequilibrium one: deviations from a description in terms of Landau quasiparticles are second 
order in U and open phase space for additional quasiparticle scattering. Moreover, the kinetic 
energy ii'kinIO = J dekekNe^,(t) already relaxes to its final value during this first episode of the 
dynamics, which is on a time scale given by ti ^ 1/ppU^. More details of the calculation can 
be found in Refs. |29l26l:i()| . 

Corrections to the second order result imply a second stage of the dynamics. Motivated by 
the explicit calculation of one (out of many) forth order term we describe them by an effective 
kinetic evolution of the nonequilibrium momentum distribution using a quantum Boltzmann 
equation. The later describes the redistribution of occupation among quasiparticle momentum 
modes due to two-particle scattering events, 

dt ^ ^fc[^'^'^(^)] - -^' E T^kp,r[N'^''it)]S';trSiek + ep- e, - e.) . (15) 

pqr 

The scattering integral Ik[N'^^{t)] contains the characteristic fermionic phase space factor 
which implements the constraints due to the Pauli exclusion principle on two-particle scattering, 

V,,,r[N'^''] = [N^^'N^^'il - iV,QP)(l - N?^) - (1 - ivQP)(l - iVQP)iVQPiVQP] . (16) 

Note that the scattering integral vanishes for arbitrary Fermi-Dirac distributions including that 
one of the zero temperature Fermi gas N'^. Hence a thermal state is an attractive fixed point 
of the Boltzmann dynamics. Since a quasiparticle picture has been established during the first 
stage we apply the quasiparticle momentum distribution of the transient state iV^^(ti) as the 
initial condition of the further quantum Boltzmann dynamics. Linearizing the scattering inte- 
gral around N^ and noting that the displacement ANf^ — N'-^^(ti) — N^ is proportional to 
U^ shows that this subsequent relaxation of the momentum distribution happens on a second 
time scale given by ^2 = l/p^W^. For small interaction strength ^2 ^ ^i- This delayed relax- 
ation behavior is reminiscent of the long nonequilibrium relaxation times in glasses. Those are 
explained by a ragged potential landscape: local minima represent transient states which are 
well-separated from the global thermodynamic ground state by energy barriers. The example 
of the Fermi liquid shows that, in many-body quantum systems, analogue bottlenecks to the 
relaxation can be imposed by particle correlations. For Fermi systems at very low tempera- 
tures this is provided by the strong quantum statistical correlations due to the Pauli principle. 
However, this restriction is not exact, and in nonequilibrium significant quasiparticle scattering 
remains effective. Therefore thermalization of the momentum distribution can be expected on 
times t > t2- The existence of the prethermalization regime in the Hubbard model and its 
subsequent relaxation has been confirmed numerically in DMFT calculations, see Sec.[7l 

Prethermalization is particularly relevant in ultracold Fermi gases: There it limits the effi- 
ciency of evaporative cooling since equilibration of the remaining atoms is very slow. However, 
successive work has shown that this delayed thermalization may be turned into a feature and 
might simplify the observation of characteristic nonequilibrium physics at zero temperature, 
for instance nonequilibrium BCS behavior |32] . 
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Fig. 3. Comparison of correlation-induced corrections to the momentum distribution function in equi- 
librium (broken line) and nonequilibrium (full line). One easily reads off the twice as large reduction of 
the quasiparticle residue at the Fermi energy in the nonequilibrium case as compared to equilibrium. 





Fig. 4. Left: Conventional scaling picture where states are integrated out around the two Fermi 
surfaces with voltage bias V (here depicted for cutoff Abg < V). Right: Flow equation approach. Here 
all scattering processes with energy transfer \AE\ < Afeq are retained in H{A{eq). 



5 Transport beyond the linear response regime through Kondo dots 

A second important class of nonequilibrium problems is posed by transport beyond the lin- 
ear response regime. While transport in the linear response regime essentially only probes the 
equilibrium ground state since the transport coefficients can be related to fluctuations in equi- 
librium, a correct description of the nonlinear response regime requires an understanding of the 
feedback of transport on the steady current-carrying state itself. In particular, many energy 
scales contribute to transport in this regime, which makes it conceptually difficult to access 
using traditional scaling techniques which focus on low-energy properties. 

An experimental and theoretical paradigm for transport beyond the linear response regime 
is realized in Coulomb blockade quantum dots in the Kondo regime, where a voltage bias V 
exceeding the Kondo energy scale Tk can easily be applied [9 lO' lllT^ . The key theoretical 
insight for understanding this situation is the observation that the shot noise of the current 
across the dot leads to decoherence, which suppresses the coherent many-particle processes 
responsible for Kondo strong coupling physics |33) . 

The flow equation method captures this physics easily since it makes the Hamiltonian suc- 
cessively band-diagonal as opposed to projecting on the low-energy subspace as in traditional 
scaling approaches. This difference is depicted in Fig. H] for transport between two leads at dif- 
ferent chemical potentials. Clearly, all current-carrying states continue to contribute in the flow 
equation scheme. It should be emphasized that there are other generalizations of conventional 
scaling approaches that also incorporate these processes |34I35I36I57] . 

Explicitly, the Kondo Hamiltonian ([S|) for two leads takes the following form 



H 



a,p,a 



(ep- 



E 



Ja'aY.^ 



'(a'p')(ap) 



(17) 
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Fig. 5. Flow equation results for the static spin susceptibility xo for nonzero voltage bias for various 
asymmetry parameters r (increasing from bottom to top: r = 1.0,1.4,1.8,2.2,2.6,3.0). The data is 
plotted as a function of the effective temperature Toff = V/{l + r){l + r^^). The dashed line is an exact 
result for the behavior in the T^s — >■ limit independent of asymmetry |39) . 



Here a', a = l,r label the two leads and the chemical potentials are given by /i;.r = zLV/2. 
The conduction band electron spin operators are defined by S(a'p')(ap) = 5 Sa '^a'p'a^ap'^apB^ 
where cr are the Pauli matrices. The couplings Ja'a describe the antiferromagnetic exchange 
interaction with the localized spin degree of freedom and are related by J/^. = JuJ^r ^-iid 
Jul Jrr = ri/Fr if thc modcl can be derived from an underlying Anderson impurity model. 

r — Fi/Fr is the asymmetry parameter of the model. 

The flow equation approach proceeds by using the canonical generator ^ and calculating 
the Hamiltonian flow ([2]) consistently including terms in third order of the coupling constant 
(details can be found in Refs. |38I39I40) ). In equilibrium this amounts to a two loop calculation 
and one recovers the well-known two loop scaling equation for the coupling constant g = p J at 
the Fermi surface 



dg 



d\nA 



fcq 



y + O,,,' 



(18) 



The calculation is actually not different in the nonequilibrium case with voltage bias V except 
that now normal ordering is done with respect to the shifted Fermi seas (see Fig. |4|). This 
changes the scaling equations significantly for Af^q ^ /Ici with 



-Trel = V 27r 



V 



1 



\niV/TK) (l + r)(l+r-i) ' 



(19) 



which is just proportional to the shot noise generated by the current. E.g. the scaling equation 
for the coupling gi at the Fermi surface of the left lead takes the form 



dgi 



dlnylfoq 



= 91 



91 



1 



1 



2 A- 



F 



rcl 



fcq 



r. 



rcl 



(20) 



One notices that if the coupling is not already too large at the scale /^ci, then the strong 
coupling flow crosses over into a weak coupling flow with g; — >■ for ylfoq — ?> 0. This permits 
the controlled evaluation of physical quantities like dynamical spin susceptibility and T-matrix 
(for details see |39I40| ). Fig. [S]for example shows the static spin susceptibility as a function 
of voltage bias, which crosses over quite accurately into the strong coupling equilibrium Bethe 
ansatz result in the limit V —>■ (see inset of Fig. [5]). The generalization to the case with 
nonvanishing magnetic field can be found in Ref. [40| , which leads to a more intricate interplay 
of different decoherence scales. 



11 
6 Nonequilibrium dynamical mean-field theory 

Dynamical mean-field theory (DMFT) 3 usually gives a reliable description of correlated 
systems in equilibrium whenever their physical properties are determined by local temporal 
fluctuations, and spatial correlations are not too important. Within DMFT, local correlation 
functions of the lattice model are obtained from an effective impurity problem which contains 
a single site of the lattice that is coupled to a bath of noninteracting degrees of freedom. This 
mapping becomes exact in the limit of infinite dimensions [2j , and it can be formulated for both 
equilibrium and nonequilibrium situations |15j . using either the imaginary-time Matsubara or 
the Keldysh formalism. Nonequilibrium DMFT has been used to study nonlinear transport 
in the Falicov-Kimball model J16I41I42I43I44) . as well as the behavior of the Falicov-Kimball 
model |45l46j and the Hubbard model [47] when the interaction is changed abruptly, or slowly, 
as a function of time. Furthermore, the method can be used for the description of time-resolved 
photoemission |48l49j and optical spectroscopy [SD] in correlated systems. 

We now briefly review the nonequilibrium DMFT formalism. A more detailed discussion of 
the technical aspects of this approach can be found in Refs. [TB] and [HT]. We then present the 
derivation of a class of closed-form self-consistency equations for the nonequilibrium case. The 
simplest self-consistency equation results for the semi-elliptic density of states, and was used 
already in several studies of interaction quenches in the Falicov-Kimball [45146) model and the 
Hubbard model [47 (cf. Sec. [71). 

The impurity problem of nonequilibrium DMFT is defined via the single-site action 



f dt Hiocit) - iY, jdt f dt' cUt)A{t,t')cAt') (21) 



ic „ Jc Jc 

on the Keldysh-contour C that runs from tmin to some time imax (i-C-, the largest time of interest) 
on the real time axis, back to imin, and finally to —il3 along the imaginary time axis [52151] . The 
first term on the right-hand side of this equation contains the dynamics due to the local Hamil- 
tonian, e.g., Hioc{t) — C/(t)n^n^ — /i(n-t--l-n^) in case of the Hubbard model with time-dependent 
interaction, e.g., as in Eq. (jl3|) . The second term describes the hybridization of the site with 
an environment that replaces the rest of the lattice, and the hybridization function A(t, t') 
must be determined self-consistently. From the action (|2T|) . local contour-ordered correlation 
functions are obtained by computing the trace {A{t)B{t') • • • ) = Tr[Tc exp{S)A{t)B{t) ■ ■ ■ ]/Z, 
where Tg is the contour-ordering operator, and real-time correlation functions can be read 
off from contour-ordered correlation functions by choosing the time-arguments appropriately. 
For the Hubbard model, the evaluation of those impurity correlation functions is the most de- 
manding part of the DMFT solution. Real-time quantum Monte Carlo methods [53] were used 
successfully 47! (cf. Sec. [7]), although the maximum accessible time is limited by the dynamical 
sign problem. On the other hand, most nonequilibrium DMFT investigations have so far been 
performed for the Falicov-Kimball model, where Monte Carlo methods are not needed because 
one can derive a closed set of equations of motion for the impurity Green functions [54| , which 
is then solved on the real time axis. 

To determine the hybridization function A{t, t') self-consistently, one must first compute the 
local self-energy S from the Dyson equation of the impurity model, 

G^[i^t+^i~A-S]-^, (22) 

where G{t,t') = —i{TcCa{t)c'l{t')) is the local Green function. Here and in the following, 
correlation functions are matrices in their contour-time variables. Matrix multiplication denotes 
a convolution along the contour C, the identity is the contour delta- function [52], and the 
operator dt denotes the time derivative. Furthermore, we restrict ourselves to the paramagnetic 
phase and omit spin indices. Because there is usually no translational invariance in time in 
nonequilibrium situations, Eq. (|22p must be solved in the time domain instead of the frequency 
domain. Next, the momentum-dependent Green function Gk = —'i{'^cCka{t)Cka.{t')} of the 
lattice model is obtained from the lattice Dyson equation 

Gk = [idt + fi - ek - S]-\ (23) 
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where spatial homogeneity has been assumed, and e^ are the band energies. The DMFT self- 
consistency is finally closed by computing the local Green function at a given site j in the lattice 
model, Gja — J2k \{j\^)\'^^ka, and requiring it to be equal to the local Green function G of the 
impurity model. If the band energies e^ are time- independent, i.e., when there is no external 
electrical field, the momentum summation can be reduced into a single energy integral 



G = Jdepie)Gie), (24) 

where p{e) = X)/c lOI^)!^*^!^ ^ ^fc) i^ ^^"^ local density of states, and G{ek) — Gk is the 
momentum-dependent Green function that depends on momentum only via the band energy 

Efe- 

The numerical solution of Eqs. (j22p and (j23p is achieved either via an explicit matrix inver- 
sion [U] (after the contour C is discretized into N time slices, and all contour Green functions, 
the derivative operator, and the identity are transformed into TV-dimensional matrices), or the 
equations are rewritten as a set of integro-differential equations of Volterra type which are 
then solved by standard numerical algorithms [5T]. Although both approaches are rather well 
behaved and can be carried on to quite long times imax (see, e.g., [46 ) the numerical effort can 
add up considerably if the k summation or e integration in Eq. (|24p requires a large number 
of fc-points |16l41j . It is therefore desirable to find cases in which Eqs. (P^ - ([M)) can be further 
simplified. Such a simplification occurs for the semi-elliptic density of states (with bandwidth 
W^A), 

V4-e2 
P(^) = —^ , (25) 

which corresponds to nearest-neighbor hopping on the Bethe lattice [55156157158) , or a particular 
kind of long-range hopping on the hypercubic lattice "59" . If the density of states ([25]) is inserted 
in Eq. (|24p . the self-energy can be eliminated from Eqs. (I^^ - (|24I) . such that one obtains a closed 
expression for the Weiss field |45) . 

A = G. (26) 

This closed form of the self-consistency reduces the DMFT self-consistency cycle to a repeated 
solution of Eq. (pS)) and the single-site problem for the local Green function G. In case of 
the interaction quench in the Falicov-Kimball model the existence of a closed self-consistency 
is crucial for the derivation of an analytical solution, thus providing the unique opportunity 
to study the limit of infinite times. We will now give a detailed derivation of this relation 
from a slightly more general statement, which in principle allows to obtain similar closed form 
self-consistency equations for densities of states other than Eq. (|25p . 

Suppose that a density of states p{e) is given and that its Hilbert transform, 

g{z)^fde-^, (27) 

J z c 

for complex frequency z satisfies the equation 

zg = l + J2 fng"- (28) 

n=l 

with an analytical function F{g) — J2'^=i fng" with real coefficients /„. Note that the coeffi- 
cients /„ can be obtained order by order by first expanding Eq. ([77)) for large z and inverting 
the resulting moment expansion into a series of z in powers of g. For example, F{g) — g^ for 
the semi-elliptic density of states (P5|) . 

Consider now square matrices Z and G that are related by 



G= dep{e)G{e) (29) 

G{e) = iZ-e)-\ (30) 
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Then, as we will show in the remainder of this section, the matrices Z and G are determined 
by the matrix equation 

oo 

ZG = 1 + ^ /„G". (31) 

which is analogous to the scalar equation ([28|. By setting Z — idt + fJ.~ I^ we see that Eqs. ((29)) 
and ([50]) correspond to Eqs. (|24p and (|23p . respectively. Equation (1^51) can be transformed into 
Eq. (P^ if we multiply it with G from the right and set 



yl-^/„+iG". (32) 

This last equation (which reduces to Eq. (0^ for F{g) = g^) provides the desired closed form of 
the self-consistency equation. We note that if F{g) is not a polynomial of finite degree, one can 
think of a suitable expansion in terms of orthogonal polynomials, e.g., Chebyshev polynomials 
|60| . Because orthogonal polynomial of G (which involve n-fold convolutions) can be computed 
recursively, this might still be more efficient that performing the e integral (j24p with one matrix 
inversion per integration point. 

To prove Eq. (PTl) for general square matrices, we multiply the equation [Z — e)G(e) = 1 
with p{e) = — Im[5(e + i0)]/7r, and integrate over e. Using Eq. (P^ . this yields 

ZG^l-- I deelm[g{>L + iQ)]G{(L). (33) 

Employing Eq. (pS)) for the scalar g, one can replace elm[(7(e + i0)] by Im[F((7(e-|-iO))], leading 
to 

ZG=\--Y,fn I delm[g{e + lO)"] G(e). (34) 

It thus remains to prove that 

G" = -- /"deIm[5(e + iO)"]G(e) (35) 

holds for any integer n > 1, which is done by induction: The initial step {n — 1) follows from 
the definition (P^)) . For the induction step, consider 

(G)"+i ^ {GYG 
(i) 1 



d€ de' lin[g{e + iO)"]Im[5(e' + iO)] G(e)G(e') 

(^) -i I & de' ImL9(e + iQT]Mg{e' + m G{e) (^7^^^^) G(e') 
(!!!) 1 [ ^ ^ , Im[.9(e + ^0)"]lm[.g(£' + ^0)] ,^, , 

^=^ -- I de G(e) (g(e + iO)"Im[g(e + iO)] + Im[g(e + iQT]g{e - zO)) 



- /"deIm[g(e + iO)"+i]G(e) 



In step (i), proposition psp is used. In (ii), the term in braces is unity, and in (iii), we use 
Eq. ([50]) . To proceed to (iv) one performs one of the two energy integrals by making use of the 
spectral representation 

TT J z- e 

The latter holds because g{z)'^ is analytic in the upper half plane. Summing up the terms in 
step (iv) and using g{e — iO) = g{e + iO)* completes the induction, and hence the proof of the 
closed self-consistency equation ([5^ . 
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Fig. 6. Fermi surface discontinuity An after quenches from U = to U < 3 (left panel) and U > 3.3 
(right panel) |47I51) . The expected thermal value for this quantity is zero, corresponding to a continuous 
momentum distribution at finite temperature. The quarter-bandwidth V is set to unity. Solid lines in 
the left panel indicate the flow equation prediction for the prethermalization plateau [29) and the 
transient behavior [25] for small values of U. 



7 Interaction quenches in the Hubbard model in DMFT 

We will now discuss the relaxation dynamics of Hubbard models after a sudden change in the 
interaction parameter from to a finite value of U, see Eq. ()13p . Since only the dynamics 
of the interacting fermions is considered, but no coupling to an external bath is present, it 
is not immediately clear whether the system will thermalize, i.e., whether after sufficiently 
long times it can be described by the thermal state that is predicted by equilibrium statistical 
mechanics. Many models of isolated many-body systems have so far been investigated (see, 
e.g., Refs. 161 62 63 64 65 66 67 68 45 29 69 70 71 72 47], and [73] for a recent review), and in 
fact thermalization is only rarely observed. In this section we will discuss results obtained with 
nonequilibrium DMFT, which show that the system can thermalize quickly for certain values 
of the interaction U |47I51| . 

For a quench from the ground state at t/ = to finite values of U the DMFT equations for 
the paramagnetic phase (Sec. [5]) were obtained in Refs. |47I51| . using the semi-elliptic density 
of states ([^5]) and real-time quantum Monte Carlo methods |S3] to obtain the Green function 
for the action (|2ip . Fig. [B] shows the jump An{t) in the momentum distribution at the Fermi 
surface as a function of time, in separate panels for small and large values of the final inter- 
action parameter U. Clearly the time evolution after an interaction quench in the Hubbard 
model depends sensitively on the parameter U. Note that An{t) remains finite for a finite time 
after the quench; in the case of a local self-energy as in DMFT this is due to the relation of 
An{t) to the retarded Green function at e = [W' . From this relation one can infer that the 
collapse of An is closely related to the decay of charge excitations that are created at the Fermi 
surface by the quench. It should be noted that as long as An{t) is finite the system is not yet 
thermalized, because a finite jump is possible in a Fermi liquid in thermal equilibrium only at 
zero temperature, whereas the quenched system has finite excitation energy. 

For quenches to U < 3 the Fermi surface discontinuity An(t) remains finite for times i < 5 
(left panel in Fig.[S]). The plateau in An{t) at intermediate times is given by 2Z — 1, where Z is 
the quasiparticle weight in equilibrium at zero temperature and interaction U (cf. Sec. 14]). On the 
other hand, the double occupation does essentially relax to its thermal value on this timescale 
|47| , confirming that the potential energy (and therefore also the kinetic energy) relax quickly, 
whereas the occupation of individual states still changes. Interestingly the prethermalization 
plateau remains well visible even for quenches to relatively large U < 2.5, even though the 
timescales V/U'^ and V^ /U^ are then no longer well separated. Not only the prethermalization 
plateaus, but also the transient behavior predicted from the fiow equation analysis |26j agrees 
well with the numerical DMFT results for U < 1.5. The prethermalization plateaus for small 
U are due to the vicinity of the integrable point a.t U — 0, as discussed further below. 
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The relaxation dynamics show different characteristics for quenclres to large U (right panel 
of Fig. [6]), namely so-called collapse-and-revival oscillations with approximate frequency 2tt/U. 
These are to due the vicinity of the atomic limit {V — 0), for which the propagator e~*^* is 
periodic with period 27r/f/-periodic [3]. For finite V these oscillations are damped and decay on 
timescales of order 1/V. For the double occupation these oscillations are not centered around 
the thermal value, but rather around a different value that can be explained using perturbation 
theory for strong coupling 47 . Somewhat analogously to the situation at small coupling the 
relaxation to the thermal state is thus delayed because the system is trapped in a metastable 
state close to an integrable point. 

Interestingly both the weak-coupling prethermalization plateau in An{t) and the strong- 
coupling oscillations disappear in a small region of interaction parameters around f/^^" w 3.2 
[57]. For quenches to values of U near C/^''" the system thermalizes very quickly. Not only 
the Fermi surface discontinuity and the double occupation relax to their thermal values, but 
in fact the retarded nonequilibrium Green function relaxes to the corresponding equilibrium 
quantity [51] . It is therefore justified to say that the system indeed thermalizes in this case, 
because a large set of observables tend to the thermal value predicted by equilibrium statisti- 
cal mechanics. Note that the statistical theory contains no adjustable parameter, because the 
effective temperature T* is determined by the energy of the system after the quench, for ex- 
ample T* — 0.84 for the quench to U = 3.3. The sharp crossover in the relaxation parameter 
is intriguing, not least because the relation to the equilibrium Mott metal-insulator transition 
is not obvious. The critical endpoint of the latter is located at Tc ~ 0.055^ |3], so that little 
of it is visible at the much higher temperatures T*. For quenches of the anisotropy param- 
eter in Heisenberg chains a remarkably similar behavior was found [72J. There the staggered 
magnetization indeed relaxes fastest to zero for quenches to the equilibrium critical value. 

It follows from these results that an isolated fermionic many-body system can thermalize 
merely under the time evolution with the interacting Hamiltonian, without requiring coupling 
to a bath. In the vicinity of an intermediate value U^-^^ the system quickly thermalizes, but 
for smaller or larger coupling thermalization has not been observed numerically for the short 
available observation times. 



8 Interaction quenches in the 1/r Hubbard chain 

For comparison we now discuss results for the one-dimensional 1/r Hubbard model, which was 
originally proposed and solved by Gebhard and Ruckenstein |74) . Its hopping amplitudes are 
given by t„ij — {—iW/2L){—l)"^~^/sin[TT{m — j)/L] with periodic boundary conditions, leading 
to a linear dispersion efc=Wfc/(27r), where W{= 1) is the bandwidth. For U > —1 this model can 
be mapped to an effective free Hamiltonian for hard-core bosons |74l75j for which the spectrum 
can be determined at once. For half-filling a Mott-Hubbard metal-insulator transition occurs at 
the interaction strength Uc ^ 1 in this model, and the Mott gap is A = U — Uc in the insulating 
phase. 

Based on this solution, the exact time evolution was obtained in Ref. |69| for a quench from 
U — (or U = oo) to a finite value of U. For the quench from U = to finite U at half-filling the 
double occupation has the noninteracting value ^ at time t = and relaxes to a new constant 
value with algebraic decay, 

m = l-'^^-'^^^ln\l^\-^^^mS^^Oa]. ,37, 



8 16U 16C/2 



2C/i2 Vt3 



In general the long-time limit of d{t) does not agree with the thermal prediction dthcrm- For 
example, for quenches to the critical value f/c = 1 the stationary value is d(t = oo) = 0.125, and 
this differs too from the thermal prediction dtherm — 0.098 that is obtained from equilibrium 
results (74) at the temperature that gives the same mean energy as the time-evolved state. 

The reason for the nonthermal steady state in this model lies in its integrability. The interact- 
ing fermionic Hamiltonian H can be mapped on an effective Hamiltonian without interactions. 
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I.e. 



iJ = ^ £„ Z„ , (38) 



where the Iq, have the eigenvalues or 1 and commute with each other (and hence with H). 
This is a rather strong case of integrabihty, because not only are there many constants of motion 
(their number is proportional to the system size), but also the energy levels can be populated 
independently (unlike, e.g., in many other models that can be solved by Bethe ansatz). As a 
consequence the fundamental postulate of statistical mechanics is not fulfilled, which states 
that all accessible microstates are assumed to be equally probable in the ensemble description 
of the equilibrium state. For an integrable system p8l) . however, all the Iq, remain at their 
values in the initial state at i = 0. That this restriction leads to nonthermal steady states has 
been observed in a variety of models [62163164165166168145] , and also in experiments with cold 
atomic gases [75] . 



9 Statistical description of nonthermal steady states with generalized 
Gibbs ensembles 

As described in the previous subsection, standard statistical mechanics cannot be expected to 
predict the steady state after a quench in an integrable system correctly. This leads to the 
question whether a suitably generalized statistical theory can make the correct prediction. The 
general procedure to construct the statistical operator in statistical mechanics is to maximize 
the entropy, and to take conserved quantities into account by fixing them on average using 
Lagrange multipliers [77]. This means that for the integrable Hamiltonian ([55)1 all constants 
of motion !„ should be fixed on average such that they yield the correct initial expectation 
value [64166] . leading to a generalized Gibbs ensemble (GGE) 

PGGE(xe-^=^°^°, (39) 

with Aq determined from (Iq)gge = (2q)o- For the l/r Hubbard chain discussed in the pre- 
vious subsection, the GGE prediction for the stationary value of the double occupation agrees 
precisely with the calculated long-time limit [SS]. However, GGEs have been found to fail for 
some observables in other models [64) . 

It is possible to formulate sufficient criteria for the steady state to be correctly described by 
the ensemble p9p . In general this depends on both the initial state \tp{t = 0)) and the observable 
A [BHl. The long-time average of ("0(^)1 A\4>{t)) is given by 

(Ay = ^(m|A|m)|(m|*'(i = 0))p, (40) 

m 

provided that the spectrum of H is nondegenerate. Note that (A) equals the long-time limit 
of (A), if the latter exists. Here the eigenbasis of the constants of motion, Zq, \m) = nia I'm), 
was used. Let us assume that the constants of motion only have eigenvalue ma = 0, 1 and 
can be represented as bosons or fermions according to Xq = aj^a^ (cf. [69] for a more general 
case)). Consider then an observable A that can be written as a linear combination of products 
of n creation operators a}^ . followed by n annihilation operators a^ . ■ The long-time average 

(A) is given by a linear combination of expectation values {Y[i=i^ai)t=o in the initial state; 
off-diagonal terms drop out due to the diagonal character of (fJOjl . On the other hand the 
GGE expectation value decouples each constant of motion from the others and instead involves 
Y[i=i{^ai)t=aj where the condition on the parameters Aq, has been used. If the expectation 
value of the product equals the product of these expectation values for each ai that occurs in 
the observable A, then the GGE is guaranteed to describe the expectation value of A in the 
steady state correctly. 
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We see that for the GGE prediction to be correct, the Xa- that occur in A must in a certain 
sense not be too correlated in the initial state, so that the expectation value of their product 
factorizes. Unfortunately the physical meaning of this condition is not obvious and certainly 
depends on the way in which the original degrees of freedom appear in the constants of motion. 
Nevertheless the GGE is certain to work for quadratic operators al^aa, and this is the reason 
why the GGE succeeds for the double occupation in the l/r Hubbard model, which can be 
expressed precisely as a sum of such quadratic operators [7S] . 

In any case we can conclude that statistical mechanics, when applied properly, does indeed 
predict the steady state of an integrable many-body system that is isolated and not coupled 
to any baths. The fact that GGEs cannot describe all expectation values correctly is actually 
not surprising. Even standard statistical mechanics fails for the expectation values of specially 
crafted operators, such as powers of the Hamiltonian or projectors onto energy eigenstates, 
which remain at their initial values and thus never thermalize. But these observables are typ- 
ically highly nonlocal and their expectation values correspond to correlation functions of very 
high order. Since they are essentially impossible to measure they are not of practical importance, 
implying no severe limitation on the applicability of equilibrium statistical mechanics. 

So far we have discussed integrable systems in the sense of Eq. (p8|) . i.e., systems of inter- 
acting particles whose Hamiltonian can be transformed into a basis in which the new effective 
degrees of freedom are noninteracting, and their occupation numbers are thus constants of mo- 
tion. However, any noninteracting Hamiltonian is of course also integrable, e.g., the Hubbard 
model with U — Q, Hq — J2k<y ^kCk^'^ka- I'^ ^'^is '-^^^ ^^^ conserved momentum occupation 
numbers c^^c^^ play the role of the conserved Tq, so that thermalization is impossible when 
quenching to Hq. 

Moreover, noninteracting Hamiltonians of course also provide useful starting points for in- 
teraction quenches, as used in Sec. |4l As discussed there, a quench to a small value of the 
interaction parameter will lead to prethermalization on an intermediate time scale [33], which 
is due the vicinity of the integrable point ai U = Q. In fact, the prethermalized expectation 
value of an observable A can be obtained by using perturbation theory and taking the long-time 
limit (provided [A, if o] = 0) [1^ . It is therefore useful to view the prethermalization plateau as 
due to the conservation of certain perturbed constants of motion Xq,, which hinder thermaliza- 
tion on intermediate time scales. It is then possible to construct a generalized Gibbs ensemble 
Pgge based on these perturbatively conserved quantities, which predicts a stationary value that 
agrees precisely with the long-time limit in second order perturbation theory, provided certain 
factorization conditions are fulfilled [78]. Namely, for an observable A = Jli^cti ^^e prethermal- 
ization plateau, which occurs for small quenches away from Hq, is correctly predicted by pgge 
if (Ili^ai) — Wi^ai)- Here the expectation value is to be taken in the (perturbed) eigenstate 
after the quench. Again it is not surprising to see that only sufficiently simple observables can 
be expected to be correctly predicted by the statistical theory. In the simplest case, the GGE 
prediction for the prethermalization plateau is correct for the observables Zq,, e.g., the mo- 
mentum distribution for quenches to small U in the Hubbard model. The agreement between 
the time-evolved state and the GGE prediction shows that prethermalization plateaus evolve 
continuously from the nonthermal steady state in integrable models, at least perturbatively. 



10 Conclusion 

For the investigation of the nonequilibrium behavior of correlated systems robust theoretical 
techniques are required, because both the interaction between particles and the time evolution 
for sufficiently long times must be described reliably. On the one hand, we discussed the flow 
equation approach and its applications to isolated many-body systems in nonequilibrium such as 
the ferromagnetic Kondo model and the Hubbard model at weak coupling, as well as nonlinear 
transport through Kondo dots. This method has the advantage that the emergence of time 
and energy scales is explicit and transparent. For example, for small interaction quenches in 
Hubbard models the method reveals that the system prethermalizes, i.e., that nonthermal 
states form on an intermediate time scale. Quantum Boltzmann dynamics then shows that 
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these nonthermal states thermahze on a longer time scale [29126) . On the other hand, we 
discussed applications of nonequilibrium dynamical-mean field theory, which maps a lattice 
system onto an effective single-site problem that can be solved numerically. The numerical 
results for interaction quenches in the Hubbard model confirm the prethermalization scenario 
and also show that thermalization can indeed occur on short timescales in an isolated many- 
body system at intermediate coupling. We compared with results for special solvable models, in 
which nonthermal steady states occur due to the presence of many constants of motion. Within 
second order perturbation theory, the statistical predictions of generalized Gibbs ensembles, 
which take these constants of motion into account, show that prethermalization in nonintegrable 
systems and nonthermal states in integrable systems can be related. 
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